﻿using System;

namespace CubicSplines
{

    /// <summary>
    /// A class for 1-d cubic spline interpolation.
    /// </summary>
    /// <remarks>
    /// The second derivatives at the node points are precalculated in constructors,
    /// if their values are not known beforehand.
    /// </remarks>
    public class CubicSpline
    {

        private double[] _X = null;
        private double[] _Y = null;
        private double[] _Y2 = null;
        private int _n = 0;
        private double _Y2_0 = 0.0;
        private double _Y2_N = 0.0;

        /// <summary>
        /// Constructor.
        /// </summary>
        /// <remarks>
        /// Makes a 'natural spline'. X[i]'s ought to be steadily increasing.
        /// </remarks>
        public CubicSpline(double[] X, double[] Y)
        {
            _n = X.GetLength(0);
            _X = X;
            _Y = Y;
            PrecomputeSecondDerivatives();
        }

        /// <summary>
        /// Constructor.
        /// </summary>
        /// <remarks>
        /// Creates a spline with known second derivatives on the
        /// boundaries. X[i]'s ought to be steadily increasing.
        /// </remarks>
        public CubicSpline(double[] X, double[] Y, double Y2_0, double Y2_N)
        {
            _n = X.GetLength(0);
            _X = X;
            _Y = Y;
            _Y2_0 = Y2_0;
            _Y2_N = Y2_N;
            PrecomputeSecondDerivatives();
        }

        /// <summary>
        /// Constructor.
        /// </summary>
        /// <remarks>
        /// Prepares a spline with values of second derivative known at the
        /// spline nodes. X[i]'s ought to be steadily increasing.
        /// </remarks>
        public CubicSpline(double[] X, double[] Y, double[] Y2)
        {
            _n = X.GetLength(0);
            _X = X;
            _Y = Y;
            _Y2 = Y2;
        }

        /// <summary>
        /// Second derivatives at the nodes.
        /// </summary>
        public double[] Y2
        {
            get
            {
                return _Y2;
            }
        }

        private void SolveTridiagonal(double[,] tri, double[] r, ref double[] u)
        {
            int i;
            double[] alpha = new double[_n];
            double[] beta = new double[_n];
            u = new double[_n];
            alpha[1] = -tri[0, 2] / tri[0, 1];
            beta[1] = r[0] / tri[0, 1];
            for (i = 1; i < _n - 1; i++)
            {
                double denom = -tri[i, 1] - tri[i, 0] * alpha[i];
                alpha[i + 1] = tri[i, 2] / denom;
                beta[i + 1] = (tri[i, 0] * beta[i] - r[i]) / denom;
            }
            u[_n - 1] = (r[_n - 1] - tri[_n - 1, 0] * beta[_n - 1]) /
                (tri[_n - 1, 1] + alpha[_n - 1] * tri[_n - 1, 0]);
            for (i = _n - 2; i >= 0; i--) u[i] = alpha[i + 1] * u[i + 1] + beta[i + 1];
        }

        private void PrecomputeSecondDerivatives()
        {
            double[,] a = new double[_n, 3];
            double[] b = new double[_n];
            a[0, 0] = 0.0;
            a[0, 1] = 1.0;
            a[0, 2] = 0.0;
            b[0] = _Y2_0;
            for (int i = 1; i < _n - 1; i++)
            {
                a[i, 0] = _X[i] - _X[i - 1];
                a[i, 1] = 2.0 * (_X[i + 1] - _X[i - 1]);
                a[i, 2] = _X[i + 1] - _X[i];
                b[i] = 3.0 * ((_Y[i + 1] - _Y[i]) / (_X[i + 1] - _X[i]) -
                    (_Y[i] - _Y[i - 1]) / (_X[i] - _X[i - 1]));
            }
            a[_n - 1, 0] = 0.0;
            a[_n - 1, 1] = 1.0;
            a[_n - 1, 2] = 0.0;
            b[_n - 1] = _Y2_N;
            SolveTridiagonal(a, b, ref _Y2);
        }

        /// <summary>
        /// Calculates the approximate value at a point.
        /// </summary>
        public double this[double xx]
        {
            get
            {
                if (xx <= _X[0])
                    return _Y[0];
                if (xx >= _X[_n - 1])
                    return _Y[_n - 1];
                if (_n == 1) return _Y[0];
                if (_n == 2) return _Y[0] + (xx - _X[0]) / (_X[1] - _X[0]) * (_Y[1] - _Y[0]);
                int i = 0;
                while (xx > _X[i]) i++;
                if (xx == _X[i]) return _Y[i];
                double h = _X[i] - _X[i - 1];
                double a = _Y[i - 1];
                double b = (_Y[i] - _Y[i - 1]) / h - h * (_Y2[i] + 2.0 * _Y2[i - 1]) / 3.0;
                double d = (_Y2[i] - _Y2[i - 1]) / (3.0 * h);
                double dx = xx - _X[i - 1];

                return ((((d * dx + _Y2[i - 1]) * dx) + b) * dx) + a;
            }
        }

    }

}
